Global stability analysis of birhythmicity in a self-sustained oscillator 
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' We analyze global stability properties of birhythmicity in a self-sustained system with random excitations. 

The model is a multi-limit cycles variation of the van der Pol oscillator introduced to analyze enzymatic substrate 
f-H ■ reactions in brain waves. We show that the two frequencies are strongly influenced by the nonlinear coefficients 

' a and /?. With a random excitation, such as a Gaussian white noise, the attractor's global stability is measured 

by the mean escape time r from one limit-cycle. An effective activation energy barrier is obtained by the slope 
of the linear part of the variation of the escape time r versus the inverse noise-intensity 1 /D. We find that the 
' trapping barriers of the two frequencies can be very different, thus leaving the system on the same attractor for 

an overwhelming time. However, we also find that the system is nearly symmetric in a narrow range of the 
parameters. 
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PACS numbers: 74.40.+k;82.20.Wt;87.10.Mn 



Some models employed to describe natural systems, such as for instance glycolysis reactions and circadian proteins 
■ rhythmics, exhibit spontaneous oscillations at two distinct frequencies. The phenomenon is known as birhythmicity, 
' and the underlying dynamical structure is characterized by the coexistence of two stable attractors, each displaying a 
^ different frequency. Being the attractors locally stable, the system would however stay at a single frequency, the one 
<^ selected by the choice of the initial conditions, unless an external source disturbs the evolution and causes a switch to the 
QQ other attractor. To investigate such process, we have focused on a particular system of biological interest, a modified van 
T— I der Pol oscillator (that displays birhythmicity), to determine the global stabiUty properties of the attractors under the 
OO influence of noise. More specificaUy, we have characterized the stabiUty of the attractors with the escape times, or the 
average time that the system requires to switch from an attractor to the other under the influence of random fluctuations. 
Such analysis reveals that the two attractors can possess very different properties, with very different relative residence 
O times. Even excluding the most asymmetric cases, the system can spend something like 10 years on one attractor for each 
second spent on the other. We conclude that although a system can be structurally biorhythmic for the contemporary 
^ presence of two locally stable attractors at two different frequencies, actual switch from one frequency to the other could 
J> be very difficult to observe. A global stabiUty analysis can therefore help to determine the region of the parameter space 
in which birhythmic behavior will be genuinely observed. 
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I. INTRODUCTION 

Self-oscillating systems exhibit limit cycles, or periodic sustained oscillations. Examples are abundant, with periods ranging 
from cardiac rhythms of seconds, glycolysis over the minutes, circadian oscillations over the 24 hours, while epidemiological 
oscillations extend even over the years uTj-yt]. Birhythmicity refers to the coexistence of two attractors characterized by two 
different amplitudes and two frequencies: depending on the initial conditions, the system can produce self-oscillations at two 
distinct periods. Such hysteretic behavior has been sometimes observed in biological systems |4]. Many more theoretical 
studies have shown the possible occurrence of birhythmicity in models of glycolytic oscillations |i5J, chemical kinetic equations 
i^, circadian proteins rhythmics and biochemical reactions ifToll . Perhaps the simplest model that exhibits birhythmicity 

is a variation of the well known van der Pol oscillator proposed by Kaiser [11] to model enzyme reactions. In such a model it 
has been shown that two attractors can coexist for some values of the parameters [11-13J, and birhythmicity is robust enough 
to enable two ifTill or more fTsIl oscillators to synchronize. The aim of this work is to adopt the Kaiser modification of the van 
der Pol oscillator llTIl as a paradigm for birhythmicity to analyze the global stability properties of the attractors under 

the influence of random excitations, i.e. the response to finite perturbations [18- 20]. In fact while local stability properties that 
refer to small perturbations of the steady state have been analyzed in Ref. flSll . global stability refers to the response to large 
random fluctuations (large enough to drive the system from one attractor to the other). Such global stability property has not 
been addressed for the model proposed in Ref. iflll - USll . and seldom investigated in birhythmic^stems (see Ref. 1I2III for an 
exception). Global stability is well studied in ac driven (and hence monorhythmical) systems 1 19l, l20ll22il23ll. for instance in 
connection with the phenomenon of stochastic resonance Ii241 or of switching between chaotic attractors ll25l l26ll . We want 
here to focus on the passage between two attractors characterized by two different frequencies, and therefore we will emphasize 
the consequences of noise driven switching on the birhythmic properties, while in periodically driven systems the frequency is 
pre-selected by the external drive. 

When noise is added, the mean time t required to escape from a basin of attraction is a useful measure of the attractor's global 
stability also for non equilibrium or oscillating systems, such as ac-driven Josephson circuits with intrinsic thermal fluctuations 
ifTsIl or with finite-spectral-linewidth ac current |27]. In the same spirit, we propose to measure the attractor's global stability 
with the mean escape time r from one stable limit-cycle attractor to another stable limit cycle attractor Escape occurs when, 
under the influence of a deterministic or random term, the system crosses the boundary of the basin of attraction (i.e. it is driven 
across the unstable limit-cycle). 

Let us remark that even if we focus on switches due to random perturbations, one could also drive the system from an attractor 
to the other by means of a deterministic or structural change. This type of switch will be not considered in the present work, 
however it is also possible from the deterministic dynamics - considering all possible paths that lead from one attractor to the 
other with the appropriate noise-dependent weight - to retrieve the escape rate ifisi I20ll28l - l30ll . 

We will show that the reason that might hamper actual observation of birhythmicity in a noisy environment is the asymmetry 
of the escape times. In such a case the system is likely to stay for a much longer time on one attractor with respect to the other, 
and therefore one would rarely observe the spontaneous transition from an attractor to the other lfl9ll20l l23^. We conclude that 
although coexistence of two stable attractors with different frequencies is a prerequisite for birhythmicity, actual observation 
might be hindered by very asymmetric stability properties of the two attractors. In other words we will consider birhythmical 
systems as bistable systems and the numerically evaluated escape times will serve as a measure of the relative stability of the 
two solutions. For a glycolytic model it has indeed been proven by means of the Fokker-Planck equation associated to the weak 
noise limit that the original system with two stable attractors (and hence with birhythmical behavior) changes structures and 
becomes monorhythmical 1,2 1,1. Our analysis arrives at a similar conclusion: the escape time from one of the attractors might 
be very large compared to the escape time of the reverse process, even by many orders of magnitude. In addition, we find that 
for some range of parameters the system is (approximately) symmetric. In this (indeed narrow) parameter space region the two 
attractors have comparable properties, and birhythmicity is more likely to be observed. 

The paper is organized as follows. In section //, we describe the self-sustained system with random excitation and the algo- 
rithm of the numerical simulations. Section ///deals with the dynamical attractors of free-noise multi-limit-cycles self-sustained 
system. We will show that birhythmicity features are not uniform in the parameter region where it appears in the modified van 
der Pol system. In section IV, we focus on numerical computed escape rates using the Box-Mueller random Gaussian generator 
algorithm |31] for numerical integration with the Euler method. The Arrhenius factor (i.e. the relation between the escape time 
T, and the noise intensity D), allows us to determine an effective activation energy barrier AUi, or the slope of the linear part of 
the variation of the escape time versus the inverse noise-intensity, as a useful method to summarize the results. The last section 
is devoted to conclusions. 
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II. THE SELF-SUSTAINED SYSTEM WITH RANDOM EXCITATION 
A. The multi-limit cycle van der Pol oscillator 

The model considered is a van der Pol-like oscillator with a nonlinear function of higher polynomial order described by the 
following nonlinear equation (overdots as usual stand for the derivative with respect to time) 

i ~ - x"^ + ax'^ - (3x'^)x + X = 0, (1) 

where a, /3 and /i are positive parameters that tune the nonlinearity. Model (1) is therefore a prototype for self-sustained systems 
and exhibits some interesting features of nonlinear dynamical systems; for instance Ref. iflil IitIi have analyzed the super- 
harmonic resonance structure and have found symmetry-breaking crisis and intermittence. The nonlinear dynamics and the 
synchronization process of two such systems have been recently investigated in Ref. lfTsllTill . while the possibility that introduc- 
ing^ an active control of chaos can be tamed for an appropriate choice of the coupling parameters has been considered in Ref. 

Eq. (|T|) describes several dynamic systems, ranging from physics to engineering and biochemistry ll33ll . In particular Eq. ([T]i 
seems to be more appropriate for some biological processes than the classical van der Pol oscillator, as shown by Kaiser in Ref. 
ll34l] . When employed to model biochemical systems, namely the enzymatic-substrate reactions, x in Eq. ([T) is proportional to 
the population of enzyme molecules in the excited polar state, the quantities a and f3 measure the degree of tendency of the 
system to a ferroelectric instability, while is a positive parameter that tunes nonlinearity li 1 3ll . 

The nonlinear self-sustained oscillator Eq. ([T]) possesses more than one stable limit-cycle solution ll34ll . a condition for the 
occurrence of birhythmicity. Birhythmic systems are of interest, for example in biology, to describe the coexistence of two 
stable oscillatory states, a situation that can be found in some enzyme reactions |35]. Another example is the explanation of 
the existence of multiple frequency and intensity windows in the reaction of biological systems when they are irradiated with 
very weak electromagnetic fields IUTI [34l [36l - [39ll . In this work we will focus on model ([T]i as a prototype for the occurrence of 
birhythmicity. 

B. The model with random excitation and algorithm for numerical simulations 

Let us consider the multi-limit-cycle van del Pol-like oscillator Eq. ([U to model coherent oscillations in biological systems, 
such as an enzymatic substrate reaction with ferroelectric behavior in brain waves models (see Ref.lfTl Ulsll for more details). 
In this case, one should include the electrical field applied to the excited enzymes, which depends for example on the external 
chemical influences (i.e., the flow of enzyme molecules through the transport phenomena). One can therefore assume that the 
external chemical influence contains a random perturbation. Therefore, adding both the chemical and the dielectric contribution, 
the activated enzymes are subject to a random excitation governed by the Langevin version of Eq. ([T]), namely: 

X - - x^ + ax'^ - I3x^)x + X ^T{t), (2) 

where r{t) is a Gaussian additive white noise ll40ll whose statistical features are completely determined by the additional prop- 
erties: 

< r(t) >= 

<T{t)T{t') >^2D5{t-t'). (3) 

The white-noise quality of F is contained in the Dirac J-function correlation The parameter D is the intensity of the Gaussian 
white noise. 

In this work we will numerically integrate Eqs. (1213b using a Box-Mueller algorithm jsTll to generate the Gaussian white noise 
from two random numbers a and b which are uniformly distributed on the unit interval [0, 1]. By introducing the new variable 



X — u, Eq. ^ can be written in the form 

X — u (4o) 
ii = fi{l - x^ + ax"^ - I3x^)u ~ X + r (46) 

The simple Euler algorithm version of the integration of equation (|4|l is given by 

Ta* = \J-^BM log{a) cos(27r6), (5a) 
x\t+At ^ X + uAt, (56) 
u\t+At — u + (/i(l — x"^ + ax* — /3x^)u — x)At + Fa*. (5c) 
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The step size used for numerical integration is generally equal to A< = 0.0001, but in some cases we have used a smaller step. 
We have also checked that averaging over as many as 200 realizations the results converge within few percents. We notice that 
there are more accurate methods to estimate the escape from a basin of attraction, or in general close to an absorbing barrier, to 
avoid the inaccuracy due to a finite sampling of the random evolution |'4l'l. However, we have carefully checked that the results 
we have obtained are independent of the step size. This has been done in two ways: halving the step size until stable results 
are reached (and with much attention to low noise intensity -0 11411]') a nd calibrating the numerical method with a potential with a 
well defined activation barrier to retrieve the Kramer escape rate 143]. 

So, although analytical treatments based on the Fokker- Planck version of the Langevin equation dU [43], the variational 
approach l|]^ri20il28l - l30ll . or faster numerical algorithms such as the stochastic version of the Runge-Kutta methods are available, 
we have preferred to use the simple procedure given by Eq. (|5]l that proved fast enough for the present project. 



S. = (a,/3) 


Analytical Amplitude 


Numerical Amplitude 


Analytical Frequency 


Numerical Frequency 


Si = (0.114; 0.003) 


Ai=2.37720 
^2=5.02638 
A3=5.46665 


Ai=2.378 
Unstable 
A3=5.464 


121 = 1.00212 

122 = 1.00113 
Sl3 = 1.0231 


Oi = 1.00015 

Unstable 

03 = 1.019575 


S2 = (0.1; 0.002) 


Ai=2.3069 
/l2=4.8472 
/l3=7.1541 


Ai=2.30265 

Unstable 
A3=7.1345 


01=0.987 
02 = 1.000113 
03=0.97123 


Oi=0.988 

Unstable 

03=0.97831 


S3 = (0.12; 0.003) 


/li=2.4269 
A2=4.2556 
/l3=6.3245 


Ai=2.4259 

Unstable 
A3=6.33918 


1=0.985 
02=0.999 
03=0.9865 


Oi=0.988 

Unstable 

03=0.988 


S4 = (0.13; 0.004) 


Ai=2A903 
A2=4.472I 
^3=5.0791 


A 1=2.48971 

Unstable 
A3=5.07739 


Oi = 1.000212 
02=1.000113 
03=0.99912 


Oi = 1.000507 

Unstable 

03=0.9989 


S5 = (0.145; 0.005) 


A 1=2.6605 
A2=3.8305 
A3=4.964 


A 1=2.65963 

Unstable 
A3=4.96336 


01 = 1.000212 

02 = 1.000113 
03 = 1.00049903 


Oi = 1.000507 
Unstable 
03 = 1.000256 


Se = (0.154; 0.006) 


Ai=2.7864 
A2=3.8821 
A3=4.2698 


Ai=2.78532 

Unstable 
A3 =4.26807 


1=0.99923 

02 = 1.000113 

03 = 1.000231 


1=0.9989 

Unstable 

03 = 1.000507 



Table 1: Comparison between analytical and numerical characteristics of the limit cycles. All data refer to the case fi = 0.1. 

In the absence of noise (F = 0), Eq. ^ reduces to the modified version of the van der Pol oscillator (see Eq. ([U), which has 
steady-state solutions that correspond to attractors in state space and depend on the parameters a, /3 and /i. Before taking up the 
subject of noise-induced transitions between dynamical attractors, we focus in the following section on the state-space structure 
of the attractors and basin boundaries in the noise-free self-sustained system. We will show that the features of birhythmicity in 
this modified van der Pol oscillator strongly depend on a and f3. 



III. DYNAMICAL ATTRACTORS AND BIRHYTHMICITY PROPERTIES 

In this Section we summarize the dynamical attractors of the modified van der Pol model ([U without Gaussian noise. The 
periodic solutions of Eq. ([T]i can be approximated by 

x{t) = Acosnt. (6) 

We recall that approximated analytic estimates of the amplitude A and the frequency have been derived in Ref. [13], and it 
has been found that the amplitude A is independent of the coefficient /i, that only enters in the frequency 57. 
It appears that, depending on the values of the parameters /3 and a, the modified van der Pol equation ^ posses one or three 
limit cycles. When three limit cycles are obtained, two of them are stable and one is unstable, a condition for birhythmicity; the 
unstable limit cycle represents the separatrix between the basins of attraction of the two stable limit cycles. We show in Fig.l 
the bifurcation lines that contour the region of existence of birhythmicity in the two parameter phase space (J3-a) {l3, 14]. The 
bifurcation line on the left denotes the passage from a single limit cycle to three limit cycles, while the right line denotes the 
reverse passage from three limit cycles to a single solution. At the conjunction, a codimension-two bifurcation, or cusp ll43ll . 
appears . The first bifurcation encountered increasing a corresponds to the saddle-node bifurcation of the outer, or larger 
amplitude cycle, while the second bifurcation occurs in correspondence of a saddle-node bifurcation of the inner, or smaller 
amplitude, cycle. The two frequencies associated to the limit cycles are very similar close to the lowest a bifurcation and clearly 
distinct at the highest a bifurcation line, as will be discuss later in detail. 



5 



Table 1 provides, for some selected sets Si of the parameters in the domain of existence of three limit-cycles on which we will 
focus our attention, the comparison between amplitudes and frequencies derived from the analytical estimate of Ref. lfTsIl and 
from numerical simulations of Eq. ([TJ. From the Table it is clear that birhythmicity is indeed present - the two stable attractors 
are characterized by different frequencies. However, the two frequencies are very similar, and in practice it might prove very 
difficult to resolve the difference. To illustrate the dynamics of the self-sustained oscillations, we report in Fig. 2 the limit cycles 
and in Fig. 3 and 4 the time dependent oscillations. In Fig. 3, the two frequencies are very similar, while in Fig. 4 we report 
the case of two clearly distinct frequencies. It is clear that for the slow oscillations (the solid line in Fig. 4, the behavior is not 
well approximated by the sinusoidal approximation (|6]l. It can also be noticed that the amplitude is still captured by the theory, 
while the agreement between the predicted and the observed frequency becomes poor at low frequencies. In fact for Fig. 4(i), 
a = 0.12, P = 0.0014, the theoretical analysis |13] predicts Ai = 2.49 and A3 = 10.89, with frequencies fli = 0.999 and 
= 0.532, respectively, in good agreement with the numerical data ili = 1.00 and Jla = 0.516. For the case of Fig. 4(ii), 
a = 0.13, (3 = 0.001, the theoretical analysis UJ gives Ai = 2.828 and A3 = 13.84, with frequencies fli = 0.998 and 
ris = 0.521, while the numerical data read fti = 1.00 and = 0.195. It is evident that the observed frequency of the large 
cycle, 0.195, is much less than the predicted value 0.521. 

In order to understand the effect of the parameters a and /3 on the dynamical states, we have simulated Eq. ([T]i to numerically 
derive the frequencies 17^; the results are shown in Table 2. For a and f3 in the white area of Fig.l, there exists only a single 
limit-cycle solution. In the gray area of Fig. 1 there are multi-limit-cycle solutions with ili / il^. Fig. 5 shows the dependence 
of the frequencies ili versus the coefficient (3 when the parameter a is fixed. In this parameter region for each value of a, the 
two limit-cycle frequencies are different at low (3 values (see Fig. 4), but converge to the same frequency when f3 increases (see 
Fig. 3). This reveals that the saddle-node bifurcation at the upper boundary of the multi-limit-cycles area in Fig. 1 occurs when 
the two frequencies are very similar Thus we conclude that birhythmicity smoothly disappears increasing /? because the two 
frequencies become undistinguishable, while the attractors are clearly distinct at the saddle-node bifurcation. 

Fig. 6 shows the dependence of 51^ versus a for different values of (3. As a increases, we move from the boundaries of the 
multi-limit-cycle area where fii = fl^ to enter the region of the map in which the two limit-cycle frequencies are different (i.e. 

So we conclude that the saddle-node bifurcation at the right hand side refers not only to the appearance of a new limit cycle, 
but also to a cycle with a definitely different frequency, and therefore in this region birhythmicity is more easily observed. In 
contrast, it is evident that it will be extremely difficult to detect birhythmicity for low a. 

IV. NUMERICAL ESTIMATE OF ESCAPE RATES AND GLOBAL STABILITY ANALYSIS 

A. Escape times from the periodic attractors 

At non zero noise intensity (D / 0), the random force causes the system to occasionally jump from one limit cycle to the 
other. The system initialized on a given limit-cycle attractor (with amplitude Ai or ^3) is forced by the random fluctuations of 
the r term in Eq. (|2]i to leave the attractor and to wander about in the neighboring state space. Escape occurs when this random 
motion drives the system across the boundary of the basin of attraction (i.e. across the unstable limit-cycle with amplitude A-z). 
The mean time r required for escape from a basin of attraction is a useful measure of the attractor's global stability. This escape 
time is analogous to the escape time of a system trapped in a minimum of the effective potential, and the escape implies that the 
random force drives the system to the other minimum of the effective potential. The activation energies shown in Fig. 7 sketch 
the escape process to be considered in the following subsection. In fact there are two metastable states: 

1. The system is trapped at the effective potential minimum in the basin of attraction of the limit-cycle amplitude Ai. Then, 
escape to the basin of attraction with limit-cycle amplitude A3 occurs when the system under Gaussian white noise crosses 
the unstable limit-cycle amplitude A2 (i.e. \x\ > A2). This can be numerically computed by choosing the initial conditions 
close to the origin. Thus, the corresponding effective energy barrier to escape from the basin of attraction with limit-cycle 
amplitude Ai to the one with amplitude A3 is called AUi. 

2. In the reverse situation, the system is trapped at the effective potential minimum in the basin of attraction of the limit- 
cycle amplitude ^3. The initial conditions are chosen outside the basin of attraction of the limit cycle Ai and far of the 
unstable limit-cycle A2. We will denote with AC/3, the effective energy barrier to escape from the basin of attraction with 
limit-cycle amplitude A3 across the unstable limit cycle with amplitude A2 (i.e. \x\ < A2) towards the limit-cycle with 
amplitude Ai . 

Fig. 7 sketches our notation and the most relevant cases: 

• Case (i): Fig. 7(i) corresponds to the case where AUi is larger than AU3. We shall see that AUi can became very large 
(depending on the coefficients a and (3); in such conditions the attractor of the limit-cycle amplitude Ai is much more 
stable than the limit-cycle amplitude A3. Thus, the system is more likely to stay on the limit-cycle attractor Ai. 
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• Case (ii): Fig. 7(ii) depicts the symmetric case AUi ~ At/3. Both attractors are equivalent and we are in a symmetric 
bistable double well. The system has approximately the same probability to stay in one or the other basins. 

• Case (iii): Fig. 7(iii) shows the case where the energy barrier AUi is less than AC/3. Here, is the reverse situation of the 
case (i), and the first attractor is less stable, the system is more likely to stay on the limit-cycle attractor A3. 

Thus, while in principle bistability occurs for all values of the parameters a and (3 in the gray area of Fig. 1, noise driven 
bistability is more likely to be observed in a narrower region of the parameter space, see case (ii). 

B. Numerical estimate of the escape rates and effective energy barriers 

Although there exists a method for the calculation of activation ener gies in non-equilibrium systems that do not admit a bona 
fide potential using the principle of minimum available noise energy |[l8l - t20l l28l - [30ll . we adopt here the indirect approach of 
computing the escape time and then we infer on the values of the activation energies. The mean escape time r is computed as 
the average over a series of trials of the time r,; required for the system to move from one attractor to the other attractor under 
the influence of noise. For each trial, integration is begun at t = with the system initialized on the attractor and proceeds by 
numerically solving the system equations with a finite difference integration method of step size At (see Eq. (|5])). The fact that 
the random motion of the system is due to a Gaussian white noise ensures that escape will occurs with probability 1 within a 
finite time 1 18]. Thus, the main question is how long the system stays in the same basin of attraction. We expect that the escape 
time is given by the inverse Kramer escape rate, or from the Arrhenius factor I.4Z1 : 

T ~ eM^UjD), (7) 

where AUi (i=l,3) is the difference between maximum and minimum values of an effective potential. 

We remark that a function plays the role of a thermodynamic potential for fluctuating dissipative systems that do not possess 
a bona fide potential llsoll if it correctly describes the asymptotic response to noise. In a sense, one reverses the Kramer logic: 
it is called effective potential a function U that gives the slope of the logarithm of the escape time vs the inverse of the noise 
intensity for low noise strength (see Eq.© ): U oc log{T/D) (for D 0). In this framework, one could regard the potential U 
as a way to summarize the behavior of the escape times. In other words it is completely equivalent either to say that the escape 
times are exponentially distributed vs the inverse of the noise (for low noise) with slope U or that the effective potential reads U. 

The relevant attractors and basins of attraction are those shown in Figs. 2. The data show that the mean escape times t 
obtained from simulations for both limit-cycle state Ai and A3, state increase exponentially with the inverse noise intensity. 
With the parameter sets Si, we find that the variation of the average escape time (on a logarithm scale) as function of the inverse 
noise intensity 1/D strongly depends to the nonlinear coefficients a and (3. For example, the sets 5*1, S4, and 5*6 correspond to 
case (i) in which the attractor of the limit-cycle ^3 is less stable than the attractor Ai. The symmetric bistable situation, case (ii) 
is observed with the set 5*5. The last case (iii) is found for the sets ^2 and S^. It is important to note that the case (ii) only occurs 
in a very narrow range, 0.08 < a < 0.09 and 0.0012 < /3 < 0.0014 lfT9ll20ll . Outside this narrow area the properties of the two 
attractors are very different. 





a = 0.07 


a = 0.08 


a = 0.09 


a = 0.1 


a = 0.12 


a = 0.13 


P = 0.004 












A (7i =0.074 
A (73 = 0.0072 


/3 = 0.003 










A(7i= 0.095 
A(73=1.656 


A (7i =0.028 
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/3 = 0.0025 
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A (7i =0.015 
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A(7i=0.25 
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A (7i =0.035 
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A(7i=0.45 
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A(7i=0.026 
A (73 = 68.2 


A (7i =0.0035 
A (73 =224 


= 0.0014 




A(7i=0.98 
A (73 =0.014 


A(7i=0.34 
A (73 =3.78 


A(7i =0.16 
A (73= 16.14 


A(7i= 0.021 
A (73= 152.3 


A (7i= 0.0017 
A(73 = 233.5 


/3 = 0.0012 




A(7i=0.62 
A (73 =2.15 


A(7i=0.291 
A(73= 11.6 


A(7i=0.13 
A(73= 17.5 


A(7i=0.104 
A (73 =308 


A (7i =0.0015 
A(73 = 791 


/3 = 0.0011 




A(7i=0.65 
A (73 =4.35 


A(7i=0.28 
A(73=27.5 


A(7i = 0.123 
A (73= 104.9 


A(7i=0.015 
A (73 = 564 


A(7i=0.003 
A (73 > 1000 


p = 0.001 


A(7i=1.3 
A (73 =0.53 


A(7i=0.52 
A(73= 10.7 


A(7i=0.25 
A (73= 16.05 


A(7i=0.11 
A (73= 105.6 


A(7i= 0.014 
A (73 >1000 


A (7i =0.0001 
A (73 > 1000 
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Table 2: Dependence of the energy barriers Af/^ in the parameters plane (a, with = 0.1. 

Fitting a straight line through the data points in the Hnear part of Eq.© and measuring its slope we obtain an estimate of 
AUi and At/3, the effective activation energies for the escape from the limit-cycle attractor Ai and A^, respectively. Since the 
effective activation energy is defined by the low-noise intensity asymptote, the accuracy of numerical simulation estimates can 
be affected if high-noise intensity points {i.e., points where the relation is not linear) are included in the fitting procedure. For 
this reason, data points for which the resulting Arrhenius factor bends have been excluded from the fitting procedure (we employ 
a test to check for linearity). Fig. 8 shows the variation of the effective energy barriers versus the coefficient /i with the set 
of parameters Si. The effective energy barriers increase when /i increases, and the behaviors strongly depend upon the set of the 
parameters Si. The scenarios mentioned in subsection IV-B can be found in the behaviors of AC/i 3 {i.e. cases (i), (ii), (iii)). The 
case (i) appears in Fig. 8 for the sets S*!, S'4, 5*6, in which the energy barrier AC/i quickly increases. Here, one concludes that 
the limit-cycle attractor Ai of the modified van der Pol oscillator is much more stable than the attractor A^ (respect to Gaussian 
white noise). The system will likely stay for a long time in the effective potential well of the limit-cycle attractor Ai, for the 
corresponding effective barrier is higher. For instance when /i = 0.5 in 5*1, we observe AC/i/ AUz — 80. The set S^ corresponds 
to the almost symmetric bistable situation, i.e. case (ii). Both effective energy barriers AC/i and AC/3 increase when /i increases 
and are comparable: the system remains for approximately the same time in the two effective potential wells. In the last scenario 
S2 and S'3, i.e. case (iii), we have a phenomenon opposed to that of the case (i): the limit-cycle attractor Aj, is much more stable 
than the attractor Ai . The system remains for a much longer time in the limit-cycle attractor Aj, because the energy barrier is 
too high, so if the noise level is large enough to cause a switch from ^3 to Ai, the same noise will drive back the system to Ay, 
in a very short time interval with very high probability. 

Let us remark that "short" and "long" might be very different lfT9ll20ll23ll . To measure the different properties, we compute 
the average persistence or residence time Pi_3 on the attractor with limit cycle amplitude Ai 3 as: 

P3 = ^^^ J = 1,3, (8) 

Ti + T3 

where ti 3 is the escape time from the first attractor Ai or third attractor A3, see Eq. (Q. For the parameters 5*1, for noise 
intensity around D = 1/20, we get P3 = 0.018, and obviously Pi — 0.982 i.e. the system will spends 1.8% of the time on the 
third attractor Aj, and 98.2% on the first attractor Ai. Decreasing the noise down lo D = 1/100, P3 decreases to P3 ~ 3.10^^. 
In other words, for any second spent on the less stable attractor A3 the system will stay for about 10 years on the most stable 
state Ai. Such a dramatic change at low noise occurs for AC/1/AC/3 ~ 50, from Table 2 it is clear that ratio between energy 
barrier can easily be much larger. 

To analyze the dynamic structure in the various areas of the chart drawn on Fig. 1, we present in Table 2 the effective energy 
barriers as a function of the coefficients a and /3 selected in the dotted rectangle of Fig. 1. When /3 is fixed and a increases, the 
effective energy barrier AC/i decreases, whereas the energy barrier AC/2 considerable increases. For example, for (3 = 0.0014, 
the effective energy barrier of the limit-cycle attractor Ai decreases from AC/i (a — 0.08) = 0.98tothe value AC/i(q: — 0.13) — 
0.0017, while the barrier AC/3 increases from AC/3(a 0.08) = 0.014 to the value AC/3(a 0.07) = 233.5. Then, there 
is a high probability that the system remains for a longer time in the limit cycle attractor A3, see Eq. (Q. A similar behavior 
is reported when a is fixed and that /3 increases. Let us note about Table 2 that for low /3 value and high a values, the case 
(iii) becomes predominant: AC/3 increases and becomes so large that we have not been able to compute such barrier even with 
simulations as long as tmax — lO^'^ normalized units. We can only estimate the barrier to be larger than 1000. 

The behavior of the effective energy barriers can be also interpreted in the following manner: the right side of the gray area 
of existence of bistable regime in Fig. 1, where the two frequencies are clearly different corresponds to the physical case where 
one of the two limit-cycle attractors, namely A3, has a very high effective activation energy while the other, namely Ai, vanishes 
because the effective potential barrier becomes zero. This process explains the passage from birhythmicity to a single limit-cycle 
attractor. 



V. CONCLUSIONS 



We have considered the characteristics of birhythmicity and the global stabiUty properties of the attractors in a self-sustained 
system. We have found that birhythmicity in a modified van der Pol oscillator is strongly influenced by the nonhnear coefficients 
a and /3: the two frequencies converge or diverge when the nonlinear coefficients are varied, leading to almost undistinguishable 
frequencies for low a and high /3. Adding a random excitation, we have found that the system crosses the boundary between 
the basins of attraction {i.e. moves across the unstable limit-cycle with amplitude A2). The mean time r to escape from one 
limit-cycle attractor to the other has been estimated in the low-noise limit, and it is proposed as a measure of the attractor's 
global stability. By considering the variation of the mean escape time r versus the inverse noise intensity 1/D, the slope of the 
linear part has enabled us to summarize the results in the form of an effective activation energy barrier which is function of the 
physical system parameters. We have found, as in other systems that exhibit noise induced switches between two attractors, that 
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the escape times can be very different lfT9il20ll23ll . so it could be difficult to observe birhythmicity for high a and low /3. We 
remark that systems lfT9ll20l 123^1 are periodically driven, and therefore monorhythmic. 

We conclude that although birhythmicity per se refers just to the occurrence of two frequencies, actual observation is subject 
to much more restrictive conditions. Our purpose is to go beyond the mere existence of birhythmicity, to show that there are 
limitations that restrict the likeliness that birhythmicity spontaneously occurs. We speculate that there might be other models 
that do possess two attractors with different frequencies, but noise driven birhythmicity is difficult to observe because of the 
different stability properties of the attractors. This might be the reason why birhythmicity has been predicted in many models, 
but rarely observed in experiments - actually there is to our knowledge just one case of clear observation of birhythmic behavior 
101 • Moreover, the switch from an attractor to another in Ref. |4j] is due to a change of the parameters, not to spontaneous 
transition from a frequency to the other We suggest that an analysis similar to that carried out in this work is therefore useful to 
ascertain the bkhythmic property in a real system. 
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FIG. 1: Parameters domain for the existence of a single limit cycle (white area) and three limit cycles 
(gray area) with = 0.1. The bifurcation line on the left denotes the saddle-node bifurcation of the outer 
or large amplitude cycle (see Fig. 2) while the right hand side contour marks the saddle-node bifurcation 
of the inner cycle. The rectangle denotes the parameter region of Table 2. 



FIG. 2: The two stables coexisting limit cycle attractors obtained by numerical integration of equation (7J 
for /X = 0.1 and the sets of parameters Si (see Table I). The thin line refers to the attractor of smaller 
amplitude (Ai) and thick line to the larger amplitude (A^). The dashed line denotes the unstable limit cycle, 
and separates the basin of attraction of the imier or smaller amplitude cycle from the basin of attraction of 
the outer or larger amplitude cycle. 
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FIG. 3: The coexistence between two regimes of noise-free self-sustained oscillations corresponding to the 
larger {solid line, A3, and frequency 0,3) and smaller {dotted line, A\, and frequency Hij limit cycles, 
the two frequencies are approximately the same. Hi ~ O3. The set of parameters is the same as for the 
attractors shown in Frg|2] see Table 1. 
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FIG. 4: Birhythmicity with oscillations at two clearly different frequencies, fli = 2fl3. (i) (a; /3) = 
(0.12; 0.0014) and (ii) (a; /9) = (0.13; 0.001). In both figures we have set /i = 0.1. 
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FIG. 5: Frequency f2j versus the parameter p for different values of a for the noise-free self-sustained 
system. The nonlinear parameter reads p, = 0.1. 




FIG. 6: Frequency SI, versus the parameter a for different values of fi for the noise-free self-sustained 
system. The nonlinear parameter reads /n = 0.1. 
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FIG. 7: Sketch of the effective activation energies AU\ and A.Uz for the free-noise self-sustained oscillator 
with multi-limit-cycles. We underline that the barrier height has clear meaning as the slope of the escape 
time\7\ while the effective potential U is qualitatively drawn only to help intuition. 
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FIG. 8: Effective activation energies versus the coefficient fi with the set of parameters Si. The thick line 
corresponds to escape from the outer cycle A3, while the dashed line refers to escape from the inner cycle, 
Ai. The parameters a and /9 are the same as in Table 1. 
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